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1. Introduction 

Gravitational Wave Astronomy has the potential to unveil the secrets of physical 
phenomena not yet understood involving strong gravitational fields. In order to 
realize this potential there are observatories/detectors either operating or being 
design/constructed that cover a significant part of the relevant gravitational- wave 
frequency spectrum. From the very low band (10~ 9 — 1CP 6 Hz), where we have the 
pulsar timing arrays to the high frequency band (1 — 10 4 Hz), where most ground- 
based detectors operate. In low frequency band (10 -4 — 1 Hz), not accessible from 
the ground, we find the future Laser Interferometer Space Antenna (LISA) |T|, an 
ESA-NASA mission consisting of three spacecrafts forming a triangular constelations 
in heliocentric orbit (see [21 IS] for details). The targets of LISA are: (i) massive black 
hole (BH) mergers; (ii) the capture and subsequent inspiral of stellar compact objects 
(SCO) into a massive BH sitting at a galactic center; (iii) galactic compact binaries; 
and (iv) stochastic backgrouns of diverse cosmological origin. 

For many of these systems it is crucial to have a priori gravitational waveform 
models with enough precision to extract the corresponding signals from the future 
LISA data stream, and also to estimate the physical parameters of the system with 
precision. Here we focus in the second type of LISA sources, namely the so-called 
extreme-mass-ratio inspirals (EMRIs). The EMRIs of interest for LISA consists of 
massive BHs with masses in the range M = 10 4 — 10 7 M Q , and SCOs with masses in 
the range m = 1 — 5OM . The SCO moves around the massive BH following highly 
relativistic motion well inside the strongest field region of the BH. The orbit is not 
exactly a geodesic of the BH because of the SCO own gravity, which influences its 
own motion, making the orbit to shrink until it plunges into the BH. During the last 
year before plunge, it has been estimated [4 that the SCO spends of the order of 10 5 
cycles (depending on the type of orbit) inside the LISA band. As a consequence the 
EMM GWs carry a detailed map of the BH geometry. This will allow, in particular, 
to test the spacetime geometry of BHs and even alternative theories of gravity (see, 
e.g. OGEE]). We can also understand better the stellar dynamics near galactic nuclei, 
populations of stellar BHs, etc. (for a review see [5]). 

All this requires very precise EMM gravitational waveforms (the precision of the 
phase should be of the order of one radian per year) . Due to the extreme mass ratio 
of these systems, = m/M ~ 10 -7 — 10~ 3 , we can describe them in the framework 
of BH perturbation theory, where the SCO is modeled as a point-like mass and the 
backreaction effects are described as the effect of a local force, the so-called self-force. 
This self-force is essentially determined by the derivatives of the metric perturbations 
(retarded field) at the particle location, which need to be regularized due to the 
singularities introduced by the particle description. The calculation of the retarded 
field has to be done numerically, either in the frequency or time domain. In this work 
we concentrate on a recent approach to self-force calculations in the time domain 
introduced in |5J [TU] for circular orbits and extended to eccentric orbits in . This is 
a multidomain approach in which the point particle is always located at a node between 
two domains, and hence it has been named the Particle without Particle (PwP) 
scheme. The main advantage is that the equations to be solved are homogeneous 
and the particle location enters via junction/matching conditions. In this way we are 
dealing with smooth solutions at each domain, which is crucial for the convergence of 
the numerical method used to implement the method, the PseudoSpectral Collocation 
(PSC) method in our case. 
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In this article we describe how to tune the numerical implementation of the PwP 
scheme to achieve high precision results with modest computational resources. We 
focus on two fronts: (i) How to change the framework introduced in jlOt ITT] to compute 
the self-force for orbits with high eccentricity, and (ii) how to pick the size of the 
computational domain in order to make our computations much more efficient, so that 
we can achieve very precise results with a modest amount of computational resources. 
A similar analysis for a different computational technique has been presented in |12j . 

The plan of the paper is the following: In Sec.|2]we describe briefly the foundations 
of self-force calculations in a simplified model based on a charged particle orbiting a 
non-rotating BH. In Sec. [3] we summarize the PwP scheme and extend the multidomain 
structure in order to allow for computations in the case of high eccentric orbits. We 
also show the convergence properties of the PSC method numerical implementation. 
In Sec. [4jwe describe how to choose the resolution depending on the mode in order to 
perform optimal calculations in the sense of computational resources. We finish with 
conclusions and a discussion in Sec. [5] where we also present some numerical results. 



2. Scalar charged particle orbiting a non- rotating BH 

The dynamics of a scalar charged particle is governed by the equation of motion of 
the associated scalar field (see, e.g. |13|): 

g^V^Or) = -4xq f dr 5 A (x, z{r)) , (1) 

J 'y 

and its own equation of motion 

dv^ dz 11 
m a ^=F^=q(g^+u^)(V u ^, «" = • ( 2 ) 

where m, q, and 7 are the particle mass, charge, and timelike worldline (parameterized 
as x^ — z^(t), being r proper time) respectively, g a " is the spacetime inverse metric 
(the BH metric, the Schwarzschild metric in our case), V„ the associated canonical 
connection, and <5 4 (x, x') is the invariant Dirac delta distribution. These two equations 
are coupled [Eq. ([TJ) is a partial differential equation whereas Eq. Q is an ordinary 
differential equation]. Eq. (JTJ) describes the retarded scalar field generated by the 
particle and Eq. Q how this field affects in turn the particle trajectory, mimicking 
the backreaction mechanism that produces the inspiral of EMPJs in the gravitational 
case (see [141 [T5] for reviews on the self- force problem) . 

Since we are dealing with a non-rotating BH, using the spherical symmetry we 
can expand the retarded field in scalar spherical harmonics 

00 l 

= £ E & m (t,r)Y lm {e,y), (3) 

t=0 m=-l 

so that each harmonic mode § em (t,r) satisfies a 1 + 1 wave-type equation: 
d 2 3 2 



dt 2 dr* 2 



V e (r)}(r& m ) = S* m 6(r-r(t)), (4) 



where r p (t) is the radial trajectory of the particle, r* = r + 2Mln(r/2M — 1) is the 
radial tortoise coordinate, M is the BH mass, and V e (r) is the Regge- Wheeler potential 
for scalar fields 

W±» + ™}, f( r ) = 1 - — , (5, 



V t (r) = f(r) 
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where S m is the coefficient of the distributional source term due to the presence of 
the particle given by 



(f.VpW) . (6) 



where tp p (t) is the azimuthal trajectory of the particle and we have assumed, without 
loss of generality that the orbit takes place at the equatorial plane 9 = ir/2. The 
singular character of this source term makes the retarded field to diverge at the 
particle location. Therefore, in order to obtain a meaningful self-force through 
Eq. (|| we must use a regularization scheme. We employ the mode sum regularization 
scheme [15, lSJ|T7J[THl[I3[2^ (see also |21[I25]). which provides an analytical expression 
for the singular part of the retarded field mode by mode. <&| m , at the particle location. 
Hence, by subtracting this singular contribution from the full retarded field modes, 
$ (which are finite at the particle location), we obtain a smooth and differentiable 
field at the particle location, <3?^ m = <3? £m — <&g™\ Adding all the harmonic modes we 
obtain the full regular field, <£> H , from which we can obtain the scalar self-force as in 
Eq. that is 



3. Formulation of the PwP scheme 



From the previous discussion it is clear that the main challenge is the computation 
of the harmonic modes of the retarded scalar field. This has to be done numerically 
and here we describe the scheme and numerical implementation that we use. We 
solve Eq. Q in the time domain. To that end we use a multi-domain framework in 
which the particle is located always at a node in the interface between two domains, 
what we call the PwP scheme. The main advantage of the PwP scheme is that it 
avoids the presence of the singularity associated with the particle in the scalar field 
equations, which is crucial for the convergence properties of the numerical method. 
The particle then enters as junction conditions of the different domains. This idea 
was already used in |23) for the computation of energy and angular momentum fluxes 
for all kinds of orbits. The scheme was used for the computation of the self-force for 
the first time in [TO] for the circular case, and extended in [TT] to the eccentric case. 
Nevertheless, the scheme as presented in [TT] can be used for moderate eccentricities, 
typically e < 0.5. The reason for this is that in the eccentric case the one-dimensional 
spatial grid needs to be changed in time in order to keep the particle at a boundary 
node. In [11] only the domains adjacent to the particle change in time, which limits 
the eccentricity of the orbits. This limitation is related to the turning points of the 
orbital motion (pericenter and apocenter), where one of the domains is effectively 
much larger than the other. In the large domain the resolution will be poor for high 
eccentricities, whereas in the small domain it will be high but this imposes a very 
restrictive Courant-Friedrichs-Lax (CFL) condition on the time step, specially in the 
case of PSC methods, where the CFL condition is proportional to 1/N 2 , being N 
the number of collocation points per domain. Therefore, pushing the scheme as it 
is for high eccentricities means to increase dramatically the computational resources 
required to obtain precise estimations of the self-force. 

Here we reformulate the PwP to allow for such high eccentricity orbits, which are 
of interest for the EMRIs that LISA is expected to observe (see, e.g. [5]). The one- 
dimensional spatial domain r* G (-co, +oo) is truncated from both sides, so that the 
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Figure 1. Structure of the one-dimensional spatial computational domain 
for a generic orbit. The trajectory is bounded to the interval between 
the pericenter (r peri ) and the apocenter (r apo ). In the setup shown there 
are three dynamical domains to the left and right of the particle, namely 

K.L.^pJ U [C+1,I.>C+1,r] U [ r m + 2,L> r m + 2,rJ ( left ) and IC+S.L > C+3,rJ U 
[ r m + 4,L> r m+4,Fj U [ r m + 5 ,L . C + S.r] ( ri § ht ) with r m+3,R = r m + 3,L = r £ • The 

figure illustrates how the dynamical domains change their coordinate size to adapt 
to the particle motion. 

global computation domain is Q — [r*,r*], where r* corresponds to the truncation in 
the direction to the BH horizon and r* corresponds to the truncation in the direction 
to spatial infinity. Now, we divide fi in a set of subdomains: 57 = uf =1 Q , where 
V a = [r* a ^ r *a, R ] with r *a,n = ''a+i.L (a = 1, . . . , d - 1). We consider two types of 
domains, static (f* L = f* R = 0) and dynamic (either f* L ^ or f* B ^ 0). For 
instance, the domains surrounding the particle, P _i = [fp_ 1L , ?"**] and J7 P = [r*,rp R ] 
(where P is the number of the domain with the particle at the left node, f2 p , and hence 
fp.j R = r* L = r*), are dynamical. Then, generalizing the setup presented in here 
we consider an arbitrary number of dynamical domains. Those domains are all around 
the particle as shown in Figure II] In the case of circular orbits (r (i) = const.) all 
domains are static. 

We discretize the spatial domain using the PSC method (see, e.g. [H]). For 
this purpose, each physical domain fl a (a = l,...,d) is associated with a spectral 
domain, [—1,1] as we use Chebyshev polynomials as the basis functions for the 
spectral expansions, which is discretized independently. Our choice of discretization 
is the Lobatto- Chebyshev grid, with collocation points given by X i = — cos(7ri / N) 
(i = 0, . . . , N). Of course, the physical and spatial domains have to be mapped and 
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this is the crucial point for our numerical implementation of the PwP scheme. For a 
given domain f2 a = [r* L , r* R l (a = 1, . . . ,d), the spacetime mapping to the spectral 
domain is taken to be linear 

X a : Tx [r* ;L) r; R ] — > Tx [-1,1] 

where T = [t ini , t end ] is the evolution time interval and 

T(t)=t, X a (t,r') = 2r '~ r ^~ rl ' R . (9) 

' a,R ' a.L 

The inverse (linear) mappings from the spectral domain to each of the subdomains Q a 
are given by 



Tx [-1,1] ^ Tx [r; L , r ; R ] 

(T,X) .— > (t,r*) 



(10) 



where 



t(T) = T , r*(T,X)\ na = ra ' R - Ta ^ x+ ra ' L Ta ' R . (11) 

The only thing left is to prescribe the time evolution of the boundary nodes of the 
dynamical domains. When there are only two dynamical domains, as in |TT], the 
evolution of their nodes is quite simple as shown above: r* ± L = const., r* R = const., 
and r* 1Ti = r* L = r*(t). Now, let us consider the situation where we have N d 
dynamical domains at each side of the particle, 2N d in total (Figure [l] shows the 
case N d = 3) . Following Figure [l] let us assume that the first dynamical domain is 
= [ r m D r m n]i which means that P = m + N d , that is 

n = K, rl, R ] U ■ • • U [C iL , r ; n J U • • • U [r* m+Ni _^, r* p ) U [r* p , r ;j 
U ■■■U[r; +JV ljL ,rJ +JV d -i, R ]U---U[r^ LI r*]. (12) 

The criterium that we use to determine the motion of the dynamical domains is to 
demand that all of them to the left of the particle have the same coordinate size (in 
terms of the radial tortoise coordinate) at any given time, and the same for those to 
the right of the particle. This completely determines the evolution of the dynamical 
domains. Then, the motion of the dynamical nodes to the left of the particle is 
described by the following equations (a = m, . . . ,m + N d — 1): 

r:, L =r:^ + (a-m) r; ~ N r ™'\ (13) 

<K=C. L + («-^+l) r; ^' L , (W) 
and the motion of the ones to the right of the particle by (a = P, . . . , P + N d — 1) 

V L =? V + (a-P) d -jjr , (15) 

Vr =r p + (a-P + l) d —^ , (16) 

where r* n L and r* +N _ l R are the fixed boundary nodes that separate the static from 
the dynamical domains. The r* -coordinate size of the dynamical domains changes 
according to the motion of the particle, as illustrated in Figure [T] 
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Figure 2. Convergence plots (log 10 |ajy| versus N) for the variable ip = r& tm . 
The upper row corresponds to a circular orbit with r = 6M (last stable orbit) 
and the lower row to an eccentric orbit with (e,p) = (0.5,8.0). The left column 
shows results for the harmonic mode (I, m) = (2, 2) and the right column for 
(i,m) = (20,0). The data have been obtained from the domain to the right of 
the particle, whose coordinate size, Ar*, is 5M for the circular case and in the 
eccentric case varies, due to the fact that the domain is dynamical, in the range 
1 — 15M. We can see how exponential convergence is achieved until machine 
roundoff is reached. 



After discretizing the domain, we need to discretize the variables and equations. 
Since nothing essential changes with respect to the description given in we 
just summarize here the main points. First we reduce the wave-type equation Q 
to a first-order symmetric hyperbolic system by introducing the following variables 
(ip,ip,ip) = (r $ , d t ip lm , d r *ij} ). Then, the restriction of these variable to every 
single domain satisfies a set of homogeneous equations, i.e. without Dirac delta 
distributions. This means that at each domain we obtain smooth solutions. This 
is the main property of the PwP scheme and of its implementation. The fact that the 
solutions at each domain are smooth preserves the exponential convergence property 
of the PSC method. At the same time, the scheme avoids the presence of singularities 
and hence, the need of using implementations that introduce artificial spatial scales 
in the problem in order to deal with the particle. We show the convergence properties 
of our computations in Figure [2j where we plot an estimation of the truncation error 
based on the value of the last spectral coefficient, a^, with respect to the number of 
collocation points, N. 

The other important point of the PwP scheme is the communication of the 
solutions from the different domains across their boundaries. This is done by using 
the junction conditions dictated by our field equations (see [251 HTJ1 ITT] ). In the case 
of boundaries that do not contain the particle this translates to impose continuity of 
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the retarded field and its first-order time and radial derivatives. When the particle 
is present the junction conditions translate into jumps in the derivatives (not in the 
retarded field itself) . The junction conditions are imposed in pratice using two different 
methods described in (i) The penalty method, and (ii) the direct communication 
of the characteristic fields. 

For the spatial discretization of our variables, U = (ip,ip,<p), the PSC method 
provides two representations: (i) the physical one, based on the values of the variables 
at the collocation points, {t/j}j = o ... jvj (ii) the spectral representation, based on the 
coefficients of the truncated expansion in Chebyshev polynomials, {a n }n=o,...,N- We 
change from one representation to the other by means of transformation matrices 
or, in a more efficient way, by means of fast Fourier transformations [SI . Finally, 
the discretization in time, the evolution algorithm, is done using a Runge-Kutta 4 
algorithm. 

4. Tuning the Domain Size to the Harmonic Number m 

In the formulation and implementation of the PwP scheme presented in [TO], [TT] it was 
already shown that precise calculations of the self-force with modest computational 
resources (less than an hour in a computer with two Quad-Core Intel Xeon processors 
at 2.8 GHz) is possible. However, as already discussed in |10l ITT] . there is a room for 
improvement in different aspects of the numerical implementation. Here, we show how 
to improve what we believe is the most important aspect: the distribution of domain 
sizes with respect to the harmonic number m. To that end, we use a phenomenological 
approach to this question, restricting ourselves to the particular case of circular orbits. 

The scalar field harmonic modes produced by a charged particle in circular motion 
are like forced oscillators, with a force term given by S lm S(r — r p ) [see Eq. Q], 
which oscillates periodically in time. More precisely, S lm oc exp(imw p f), where 

Lu p = M x l 2 r v 3 ^ 2 is the particle's angular velocity. The potential V t decays quite fast 
away from the location of its maximum, so that the field looks like a monocromatic 
wave in the regions where the potential is not important. The wavelength of these 
waves is A m = Xi/m, where Ai = 2tt/lu p is the wavelength of the m = 1 modes. 
Therefore in a domain f2 OJ with coordinate size Ar* = r* R — r* L , modes with different 
harmonic number m require different resolutions, understood as the ratio Ar* /N, in 
order to resolved them with the same accuracy. For instance, given two modes $> im 
and $ £ m such that m < m' , i.e. with A m > \' m , in a given domain fl a we will find 
more wavelengths of the mode <E> £ m than the mode 4> fm . This means we need more 
resolution for the mode <E> £ m than for the mode <& lm in order to achieve comparable 
accuracy. 

Here, we can change Ar* and the number of domains d (if we reduce Ar* we 
need more domains and the converse), or TV, or both. The question here is what is the 
optimum way of setting the resolution for each m-mode, in the sense of minimizing 
the computational resources. Increasing the resolutions means to increase the number 
of computer operations, and in this sense it is not the same to increase the number 
of domains, decreasing their size (see [10 for a discussion of this point), than to 
increase the number of collocation points (using FFT instead of matrices, the number 
of operations in the relevant part of the algorithm scales with N as N In N) . In order 
to investigate this question we have conducted a series of simulations to study the 
behavior of different modes of the retarded field $ fm under changes of the coordinate 
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domain size Ar* and the number of collocation points N. The results for values at 
the particle location for the modes $ 8 ' 2 and $ 8 ' 8 are shown in Figure H As we can 
see, it is crucial to choose carefully Ar* for a fixed value of N, otherwise we can 
find ourselves in one of the following situations: Either Ar* is too big and then, N 
collocation points are not enough to resolve all the modes, or Ar* is too small and 
we are using too many collocation points, i.e. we are performing too many numerical 
computations to resolve the modes. 




5 10 15 20 25 30 5 10 15 20 25 30 

N N 

Figure 3. This figure shows values of the <I> 8 ' 2 (left) and cE> 8 ' 8 (right) modes of 
the retarded field at the particle location as computed from the domain to the 
right of it (Hp), with respect to the number of collocation points N. Each line 
correspond to a different coordinate size of the domain where the calculations are 
done: Ar*/M = 1,5,10,20,40. It shows how the values of the modes converge 
as Ar* is decreased and N is increased. It is remarkable to realize that for small 
coordinate size the convergence is reached with few collocation points. 

From our simulations we have found that adjusting the size of the domains to 
the wavelength of the modes of the retarded field, i.e. Ar* ~ A m , the values of 
the field converge and are resolved for a minimum of collocation points given by 
N ~ Ar*/(2Af) + 25. This provides us with a rule of thumb for the choice of the 
domain size and number of collocation points as a function of to: Ar* (to) ~ X m and 
N(m) = [Ar*/(2M) + 25] (where here [x] denotes the integer closest to x). 

Another strategy that has been used in the literature (see |25l IT2"] ) is to truncate 
the mode series for the self- force at a given I — £ T and then, to estimate the remaining 
infinite number of terms by using the known large- £ series The idea is to fit 

the coefficients of this series by least-squares fitting to a subset of the numerically 
estimated self-force modes (i.e. with £ < £ T ). We have not used this procedure in this 
paper because our aim is to tune our numerical scheme, but it can certainly be used 
in order to complement and improve the results that we obtain. 

5. Conclusions and Discussion 

The PwP scheme introduced and developed in [51 QUI E] provides a framework for 
precise and efficient computations of the self-force in the time domain. In this paper 
we have investigated how to tune the method in order to increase its efficiency and 
accuracy. We have also extended the PwP scheme in order to make it suitable for 
computations of the self- force on highly eccentric orbits. 

We have studied how to distribute the domain sizes and the number of collocation 
points so that we allocate the optimum resolution to each harmonic mode. This is very 



Tuning Time-Domain Pseudospectral Computations of the Self-Force 



10 



important as the resolution requirement increases with m, despite the fact that high 
^-modes contribute less to the self-force than low I ones. This means that in order to 
improve the accuracy of the self-force we must have a good control of the resolution, 
otherwise modes with high m (and hence with high t) will limit the precision despite 
not contributing much to the self-force. 

Using this information, we have conducted a series of computations of the self- 
force for the case of a scalar charged particle in circular (geodesic) orbits around a non- 
rotating BH. The computational parameter space that we have explored is described 
as follows: The total number of domains that we have used is in the range d = 20 — 43. 
The coordinate size of the large domains (the ones far from the particle and near the 
boundaries) is in the range Ar* = 50 — lOOAf . The number of collocation points has 
been fixed to N = 50. The largest £, £ max , considered in the computations is in the 
range 4, ax = 20 - 40. 

Adapting the resolution for each m-mode also allows us to adapt the time step as 
the CFL condition is proportional to the minimum coordinate physical distance (as 
measured in terms of the tortoise radial coordinate) between grid points (and inversely 
proportional to the square of the number of collocation points). This means that the 
number of time steps required, N t , for a simulation is going to be proportional to m 
as: N™ — m A™ =1 . This is assuming that we use the same number of domains for 
all modes, but it is clear that modes with low m would need less domains than modes 
with high m, and therefore, this is another source of reduction of computational time. 

We have been able to obtain very precise values of the self-force for a wide range 
of values within the parameter space. For instance, in the case of the radial component 
of the regular field, , the only one that for the circular case needs regularization, 
we have obtained values like <f>^ = 1.677282 x 10~ 4 q/M 2 , which have a relative error 
of the order of 5 • 10 -5 % with respect to the values obtained in [25] using frequency- 
domain methods. The relative error is always in the range 5 • 10 -5 % — 5 • 10~ 3 %. This 
constitutes a significant improvement with respect to our own previous estimations 
presented in [ID], where the relative errors quoted for $^ were of the order of 0.1%. 
The typical time for a full self-force calculation in a computer with two Quad-Core 
Intel Xeon processors at 2.27 GHz is in the range 10 — 15 minutes, which is a very 
significant reduction with respect to our computations in |10| . specially taking into 
account that we have also improved the precision of the self- force computations. 

The calculations we have presented can be further improved in terms of 
computational time, and perhaps in accuracy by exploring techniques to bring the 
boundaries closer to the particle without degrading the accuracy of the field values 
near it. This can be done either by improving the outgoing boundary conditions 
(see, e.g. [26J or by using some sort of compactification of the physical domain (see, 
e.g. [27]). Another possibility for making the computations faster (although this does 
not decrease the CPU time) is to parallelize the code and use computers with many 
cores. This is in principle a simple task as the different modes are not coupled. In any 
case, the next step in this line of work is to perform a similar phenomenological study 
for the eccentric case, where not only the azimuthal frequency is important, but also 
the radial one, which is absent in the circular case. 
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